clear all;
close all;
clc;

% create time vector
maxt = 20;
deltaT = 0.05;
t = 0:deltaT:maxt;

% create space vector
L = 10;
deltaX = 0.05;
x = 0:deltaX:L;

% create vector of initial conditions (sine wave)
% and set as time steps -1 and 0.
% solutions for every time step will be stored in the columns "us"
maxh = 2.5;
f = transpose(maxh*sin(pi/L * x));
us = zeros(length(x),length(t)+1);
us(:,1) = f;
us(:,2) = f;

% create tridiagonal matrix by first creating 3 vectors containing diagonal
% matrix values then using spdiags command or a sum of diag commands
B = zeros(length(x),3);
B(:,1) = ones(length(x),1);
B(:,2) = -2*ones(length(x),1);
B(:,3) = ones(length(x),1);
d = [-1,0,1];
A = spdiags(B,d,length(x),length(x));
% A = diag(B(1:(m-1),1),-1) + diag(B(:,2),0) + diag(B(1:(m-1),3),1);


% start clock to measure computation time.
tic();

% compute solution for each time step.
for ii = 2:1:length(t)-1
    us(:,ii+1) = 2*us(:,ii) - us(:,ii-1) + deltaT^2 / deltaX^2 * A * us(:,ii);
end %for

% stop clock and output computational time
toc();

% setup figure parameters to overwrite last plot at every plot command
fig1=figure(1);
set(fig1,'NextPlot','replacechildren')

% turn on pausing to allow human visualization of "movie"
pause on;

% plot solutions at each time step.  The axis command is used to maintain 
% plot window dimensions.
display('creating movie');
for ii=2:length(t)
    plot(x, us(:,ii));   
    axis([0 L -(maxh+0.5) (maxh+0.5)]);
    pause(0.05);
end


